Effective algorithm of analysis of integrability via the Ziglin's method. 
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In this paper we continue the description of the possibilities to use numerical simulations for 
mathematically rigorous computer assisted analysis of integrability of dynamical systems. We sketch 
some of the algebraic methods of studying the integrability and present a constructive algorithm 
issued from the Ziglin's approach. We provide some examples of successful applications of the 
constructed algorithm to physical systems. 
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I. INTRODUCTION 

In this paper we describe the continuation of the research aimed at application of numerical methods to analysis 
of integrability of dynamical systems. In the first paper of the series ([l|) we have described the way of revealing the 
obstructions to real integrability in the Liouville- Arnold sense (0) by analyzing the topology of the phase space of 
the system. We have also discussed the possibilities of extending the approach to parametrized families of dynamical 
systems with the goal of searching for the regions of possible integrability. 

The current paper is devoted to analysis of algebraic properties of integrable systems. We mention some recent 
results related to the study of complexified systems and systems of variational equations in particular. We pay special 
attention to the results of S.L. Ziglin ([3() on the meromorphic integrability. We point out the major difficulties of 
application of these results to the study of integrability and propose an effective algorithm of a computer assisted 
construction significantly extending the range of their applicability. We give some examples of application of the 
method to dynamical systems having physical origin, and also mention some possible purely mathematical outcome. 



II. ALGEBRAIC OBSTRUCTIONS TO INTEGRABILITY 

Since for a given dynamical system there is no general approach for studying the existence of the sufficient number 
of arbitrary first integrals, a natural idea to restrict the class of first integrals comes out. A rather detailed review of 
the appropriate methods can be found in some more recent methods are also well explained in Here we will 
only sketch some of them, that are important to understand our motivations and the results presented in this paper. 

One can probably say, that the first algebraic method for analysis of integrability is the approach of S. Kovalevskaya 
developed in Q by H. Yoshida for studying the polynomial first integrals of dynamical systems. Let us note that 
this restriction is rather reasonable since for mechanical systems with polynomial hamiltonian functions the natural 
integrals (energy, angular momentum etc.) are of this class. H. Yoshida applies this method to study the Euler 
equations as well as some hamiltonian systems having physical origin. 
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A. Variational equations, monodromy group. 

Some more recent methods of analysis of integrability deal with the complexified systems of differential equations. 
The key observation is that rather often for classical completely integrable systems the first integrals can be continued 
to the complex domain of the canonical variables depending on complex time remaining in the class of holomorphic 
or meromorphic functions. And branching of the solutions of the hamiltonian equations can create an obstruction to 
the existence of such first integrals. This idea permits in particular to find the relations between the parameters of 
a system necessary for complete integrability, that is select more regular systems from the family of a priori similar 
ones. 

Following partially @ let us formalize the above statement. The following approach actually dates back to A. Lya- 
punov who proposed to study the system of variational equations, it was developed by S. Ziglin in [3j who revealed 
the relation of the structure of the monodromy group to integrability. 

Consider in C" (or any n-dimensional complex manifold) a system 

x = v(x), (1) 

where the dot denotes the derivative with respect to the complex time. Denote xo(t) - a particular solution of ([T]) 
and consider a small perturbation x = xo + £ of it. Note that here and in what follows the variables like x, v will be 
vector- valued (in C", in a complex manifold or in the appropriate tangent bundle), and those like £ can be identified 
to a respective tangent vector field in the whole phase space. Plugging x into ((T|) one obtains 

t = A(t)t + ..., A=j£( XQ (t)). 

Then the linearized system 

i = A(xa(t))£ (2) 

is called the system of variational equations along the particular solution xo(-)- We write A(xo(t)) instead of just 
A(t), to stress the fact that the matrix of this linear system depends on time only implicitly via its dependence on 
the particular solution. If the system @ has multivalued solutions then the system (JTJ also does. We will be however 
interested in a more subtle property of these systems, namely branching of the solutions of the system of variational 
equations ^) when moving along the Riemann surface of the particular solution xn of ([Tj). 

Thus we will consider a system of n linear equations ([2]) , where the entries of the matrix A are holomorphic functions 
defined in a connected neighbourhood of the Riemann surface xq. Locally for any given value of £(io) — Co there 
exists a unique holomorphic continuation of the solution ^ . One can continue it along any path in xq but the result 
in the generic case need not be single- valued. Let 7 be a loop (closed path) starting from x(to) parametrized by some 
path on the complex plane. Constructing the continuation of the solution £(i) of |2} defined in the neighbourhood 
of x(to) along 7 we obtain another solution £*(t) of © defined in the same neighbourhood. Since the system @ is 
linear, there exists a complex n x n matrix T 7 such that £*(f) = T 7 £(t) for any £. Branching of the solutions of © 
corresponds to T 7 ^ id. 

The set of all matrices M. = {T 7 }, corresponding to all the loops 7 in xq, forms a group called the monodromy 
group of the linear system ([2]) along the solution xq. The matrices certainly depend on the choice of the base point 
on Xo but all the groups M. are isomorphic. It is important to note that two nomotopic loops produce the same 
monodromy matrices, therefore it is sufficient to consider only one representative 7 from each homotopy class. 



B. Results of S. Ziglin 

Let f(x) be a first integral of (fTJ. One can write down the Taylor series for J(xq + £) 

/(Z0+0= X^mUo (£>*), 
rn>0 
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where F m is the homogeneous form of degree m, single- valued on the Riemann surface xq- It is clear that a non-zero 
homogeneous form F m of the lowest degree (m > 1) is the first integral of the system of variational equations ©, and 
thus invariant under the action of the monodromy group. 

F m (n,to)=F m (Z,t ), TeM. 

If f(x) is meromorphic, it can be represented as the ratio of two holomorphic functions: q ffi . Then the analogous 

P (z) 

invariant is the ratio of the lowest degree forms corresponding to P and Q: q1^ z ) ■ This is a strong condition on 
invariant functions, but even more important, the existence of such invariants imposes serious restrictions on the 
structure of the monodromy group. 

Without going much into details let us just mention that for what follows it is important to exclude the cases 
when all the matrices of the monodromy group have trivial eigenvalues. This happens when the system possesses 
so-called symmetry fields; if this is the case one needs to perform the reduction of the system to normal variational 
equations. Let us also note that for the case of hamiltonian systems, which is of particular interest for us in the 
context of analysis of integrability, the elements of the monodromy group define symplectic affine transformations. 
In this case the eigenvalues of the monodromy matrices split into couples of mutually inverse complex numbers. The 

transformation is called non-resonant if its eigenvalues (Ai, A^ 1 , . . . , A p , A" 1 ) satisfy the relation A* 1 Xp p = 1 if 

and only if all the ki vanish. 

Having defined the notion of the monodromy group we can state the following theorems: 

Theorem 1. (Ziglin's lemma) Let the monodromy group of the curve xo contain a non-resonant transformation 
g. Then the number of meromorphic first integrals of the hamiltonian equations in the connected neighbourhood of Xq 
functionally independent from H is not bigger than the number of rational invariants of the monodromy group. 

Theorem 2. (S. Ziglin) Let the monodromy group Ai of the curve Xq contain a non-resonant transformation g. If 
the hamiltonian equations admit (n—l) meromorphic first integrals in the connected neighbourhood of Xq functionally 
independent from H then any other transformation g' € Ai preserves the fixed point of g and maps its eigendirections 
to eigendirections. If moreover no set of eigenvalues of g' forms a regular polygon on the complex plane (with the 
number of vertices >2), <?' preserves the eigendirections of g (i.e. commutes with g). 



C. Differential Galois theory 

Despite some important examples, the analytical computation of the monodromy group is a complicated task, solved 
only for some particular systems. In this context it is worth noting that there is an alternative approach proposed 
by J. Morales-Ruiz and J. -P. Ramis (@, H|) following the scheme of S. Ziglin. It establishes the relation between the 
properties of the differential Galois group and integrability. Not defining this group here (for the definition see e.g. 
9]) let us only state the main result. 

Theorem 3. (Morales-Ramis) Let the hamiltonian system be completely integrable, then the identity component 
of the differential Galois group of the system of normal variational equations along any particular solution is abelian. 

When restricted to the case of complete integrability, i.e. the existence of a sufficient number of independent first 
integrals in involution, the theorem of Ziglin is the consequence of the one of Morales Ramis. The latter is considered 
to be stronger also in the sense that it permits to show non-integrability in some cases when the Ziglin's approach 
did not give any answer. The explanation of this fact is that for a fixed particular solution of ([T]) the differential 
Galois group contains the monodromy group, and therefore has more sources of non-commutativity. In this paper we 
however restrict our attention to the Ziglin's approach, leaving the results of Morales-Ramis to further studies. 



III. EFFECTIVE ALGORITHM OF APPLICATION OF THE ZIGLIN'S METHOD. 



We can formulate the standard way to analyze the integrability which is essentially based on the results reviewed 
in the previous section. The major steps are the following: 
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1. For a given complexified system, of differential equations ([T]) construct explicitly a solution xq, which as a 
function of complex time is viewed as a Riemann surface. 

2. Write down the system of variational equations along xo- Perform a reduction of this system to the system of 
normal variational equations using known first integrals or symmetry fields. 

3. Localize the singularities of the particular solution x . 

4. Construct the monodromy matrices, obtained by going along the loops around the singularities obtained in 3. 
- they generate the monodromy group. 

5. Make a conclusion about the presence of the obstruction to integrability, based on the commutation relations 
between the matrices obtained in 4. 

Let us comment on some subtleties of application of this method. The method is obviously not intended to prove 
integrability, and there are examples ([HHI) when the monodromy group is trivial, but the system is still non-integrable. 
But it can be also used to single out the relations between the parameters of the system, when integrability is possible. 
It is usually done, when the variational equations reduce to some "classical" well-studied systems. In principal, one 
can argue for integrability, when the monodromy group is commutative for any particular solution, but this is not 
worth the efforts, since knowing explicitly all the solutions of the initial system ([1]) one can perform a more detailed 
qualitative study of it. So the method can be developed mainly to search for the obstructions to integrability. 

We can easily see two major difficulties of application of this method. First, one needs to construct an explicit 
particular solution of a system of differential equations. And second, this solution should be on one hand simple 
enough, so that the the monodromy group could be computed, and on the other hand, it should be non-trivial, so 
that the computed group could contain sufficient number of sources of non-commutativity. We employ numerical 
methods to overcome these difficulties, more precisely the idea is to compute numerically the generators of the 
monodromy group along a particular solution which is also obtained numerically. 

For a given value of complex time t the value of the solution x(t) can be computed numerically without any difficulty. 
One needs just to remember, that the obtained particular solution x(t) should be considered not as a function of a 
point t G C in a complex plane, but as a function of (the homotopy class of) a path going to t from the initial point to. 
We can always write down the right hand sides of the variational equations £ = A(x(t))£ with the matrix A depending 
on an arbitrary particular solution x. The difficulty is that if we don't have the analytical expression of x[t) we can 
not plug it into the variational equations, but with the above remark it does not make much sense in the general case. 
That is why we use a natural approach of solving the variational equations in parallel to the initial system along any 
path that interests us. 

So, using the observations about the structure of the monodromy group from the previous section l|II A[) let us solve 
numerically the following system 

x = v(x) 

3 = A(x)E, (3) 

where the first line is identical to the initial complexified system of n differential equations (JTj) , the second one is the 
matrix equation with A being the matrix of the system of variational equations depending explicitly on x, and 3 - 
unknown n x n matrix. For initial data we take Xo(to) - an arbitrary point in the phase space and 3(io) = id n - unit 
matrix. Going around the loop on the solution x we obtain in 5 precisely the monodromy matrix corresponding to 
this loop. Thus going around all the singularities of the solution x we can construct the whole monodromy group. 

It is important to understand that these loops are paths on the Riemann surface of the solution and not just on 
the complex plane. This is the major difficulty that one faces when applying the procedure: given the form of the 
equations © we are forced to integrate them against the complex time, but the object of interest for us is the behavior 
of the solution x. More precisely, when we go along a path in the complex plane and neither x nor 5 have non-zero 
variation, this path does not produce a non-trivial monodromy generator. If 3 has non-zero variation, one needs to 
check that the corresponding values of x[t) have returned to the initial value. Only in this case 3 is the generator 
of the monodromy group. If x(t) has not returned to its initial value one needs to continue following the path. This 
difficulty is related to the fact that in general the solution x can not be parametrized by the complex time t, i.e. going 
around the loop in C does not always produce a closed loop on x and the parametrization used in the original work of 
S. Ziglin should be considered as the parametrization of the loop on x. It is also clear that if x does not return to its 
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initial value after a finite number of going around a loop in C (that is the topology of the singularity is logarithmic), 
such a loop does not correspond to any matrix in the monodromy group. But such an infinite branching by itself 
can be an obstruction to, say, analytic integrability, at least if the corresponding energy level is compact in the phase 
space. If all the branching points of x are of finite order (such a system satisfies the generalized Painleve property), 
then the outcome of the procedure is the set of monodromy group generators. Having that, to prove non-integrability 
it is enough to find a couple of non-commuting matrices among them. 

There is a couple of issues also worth being commented on. First, there is a difficulty which is more technical than 
conceptual, it is related to the localization of the singularities of the solution x (step 3 in the method above). Certainly 
it is impossible to go around all the loops in C, so we have to restrict the analysis to some compact domain and go 
around the points of some finite grid. That is we do not try to construct the whole monodromy group, but only its 
subgroup, which is however usually enough to reveal the obstruction to integrability. Second, we have not discussed 
here the issue of resonant transformations. But, given the symplectic nature of the monodromy transformations this 
problem arises only from relatively large size of the system ((I) in question. In this setting let us also note that one 
does not have to perform the reduction (in step 2 of the method) of the system of variational equations , it is only 
necessary to know that it can be performed. It is also important to mention that since the final action (step 5) is the 
verification an open condition of the commutator non- vanishing, it is perfectly correct from the point of view of the 
accuracy of numerical integration. 

Summing up, let us formulate the effective algorithm of analysis of integrability via the Ziglin's method. 

1. Write down (analytically) the system of equations ^ not fixing the particular solution. 

2. Choose a bounded domain in C and a finite grid of points in it with a distinguished point to- 

3. For each point choose a loop going around only it and starting at t^. Integrate numerically the system along 
this loop taking S(to) = id as the initial conditions. Three cases are possible: 

i. x and 5 returned to initial values - this point gives a trivial transformation from the monodromy group. 

ii. The value of x did not return to the initial values (within a given precision) - continue integrating around 
this loop. If x does not return to the initial values after a sufficiently large number of loops, one needs to 
analyze the density of the trajectory in the phase space. 

iii. The values of x returned to the initial values after a finite number of loops, but of S did not - store the 
matrix S, it is one of the generators of the monodromy group. 

4. Compute the pairwise commutators of all the matrices obtained in 3. iii. If there are non- vanishing commutators 
make a conclusion about non-integrability; if not choose another initial value of x(to) in step 3. 



IV. APPLICATION 



In this section we apply the developed algorithm to some systems having mechanical interest. We start with the 
example that served one of the motivations to favor the approach of Morales-Ramis in comparison to the Ziglin's one 
- the system described by the Henon-Heiles hamiltonian. Within the framework of our algorithm we can also study 
the systems we were interested in while describing another approach in namely the pendulum- type systems and 
satellite dynamics. 



A. Henon-Heiles system 

The Henon-Heiles hamiltonian describes a very simple model of a star moving close to the galactic center. It reads 

Hh = \(p\+pl) - ql{A + qi) - ^qf. 

In Q (referring also to original works @, HI) the monodromy group and the differential Galois group were constructed 
for the hamiltonian equations governed by Hh using a rather simple particular solution satisfying q2 = P2 = 0. 
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For this solution the monodromy group is trivial and does not obstruct integrability, while for A = the theory of 
Morales-Ramis permits to show that the system is non-integrable. But for A =^ and A ^ one can only state 
non-integrability for j ^ k ( k + 1 ) ^ k *E Z. For A = the approach does not give any result. Later Morales and Ramis 
performed a more detailed study, showing that the question of integrability is open only for A = 1, 2, 6, 16. Using the 
fact that one is not forced to be restricted to the above particular solution in our numerical approach, we can study 
the remaining cases by our method using a more complicated one. 

For example for A = 1, A = 0.25 consider the initial data (qi,q2,pi,p2) = (1,-0.4,-1.25,-0.3), to = 1 (the 
figures are chosen arbitrarily not to obtain a vanishing solution for q2,P2)- The commutator of matrices obtained by 
going along a loop around the points (0.2 + 2.5*) and (0.2 — 2.5*) is equal to 

0.390306 - 0.912711i 0.657898 - 3.4626* -1.09917 + 2. 5655? 0.65533- 2.22299* 

-0.936539+ 2.30099* -0.636668 + 2.59892* 1.31112 - 1.86606z -1.72067 + 7.99614* 

0.314644- 0.7603* 0.273818 - 1.21218* -0.51959 + 0.878226* 0.575715 - 2. 514* 

0.463812- 1.0788* 0.797168- 4.19885* -1.33512 + 3.13256* 0.765952- 2.56444* 

that results in non-integrability of the system. 



B. Triple pendulum 



A triple pendulum is a system of three mass-points (described by the radius- vectors r^) connected by weightless 
inextensible rods. We consider a free planar motion of this system which within the Lagrangian formalism can be 
described by a system with constraints of the form 

r? - J? = 0, 
(r< - iw) 2 - l\ = 0, i = 2,3. 

We are not going to describe the formalism in full details here since we have already sketched it together with the 
motivations in [lj (cf. also references therein). Let us only recall that using a convenient parametrization by angles 
01,02 between the segments of the pendulum the system can be reduced by Routh transform to two degrees of 
freedom. One technical difficulty is that we need to have an explicit form of the system of variational equations, 
which is in this case much more complicated than in the previous example. Luckily, we can use algorithms of symbolic 
computation to obtain it; in this case we used the Sage ([l3|) software package. 

Turning to the results, for the initial data (/3i, $i, 02, $2) = (0.3,-1,-0.15,0.5), to = 1, the commutator of the 
matrices obtained by going along the loop around the points (0.2 + 0.5z) and (0.2 — 0.5*) (each of them six times), 
reads 

/ 0.6170- 0.6190* 0.8124 + 0.8314i -0.2373- 0.3366* -0.0180 + 0.3343* \ 
-0.0658- 0.1425* 1.0420 + 0.0975* -0.0447+ 0.1958* -0.0189- 0.3878* | 
0.0206 - 1.0089* 0.2544+ 1.5232* 0.1464 + 0.5693* 0.8635- 1.4386* ' 

\ -0.0047 + 0.0235* -0.0163- 0.0706* -0.0227- 0.0896* 1.0258 + 0.1094* / 

That is the system is meromorphically non-integrable which is in perfect agreement with the results of [l|. 



C. Satellite dynamics 



Let us now consider another example already mentioned in [l[ - the motion of a dynamically symmetric satellite 
along a circular orbit f[ll|). Using again the Routh reduction procedure one can describe the dynamics by the 
hamiltonian 

2 sin * £ 
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Consider the trajectory starting at ip = 0, 9 = = 0.1, pe = 0,t = and continue it along the loop around the 

points t — 4.8 + 0.8i and t = 4.8 — 0.8i (two times around each of them). The respective commutator reads 

/ 8849.8+ 13. 2915? 37.9467 - 126. 07H 2044.45+35.4031? -1843.31 - 125.866i\ 

-9456.28- 239.498i 311.834-62.6925? -2350.67 - 37.0663i 1972.34 + 197.277z 1 

-34540.9 - 527.522* 596.362 - 53.1627i -8340.91 - 82. 1358i 7205.45 + 624.613i ' 

\ 4032.64- 556.427i 1615.13 - 971.49i 177.875 + 135. 104i -820.725+ 131.537?/ 

That is the system is also meromorphically non-integrable. An interesting observation here is that the configuration 
of the system is parametrized periodically by the angles ip, 0, that is we didn't have to make them return exactly 
to the initial condition, but only modulo 2ir. Another interesting feature of the system revealed by the numerical 
experiment, is that the complexified dynamics of it is much more sophisticated than the real one. That is if the 
system is indeed locally integrable as we have conjectured in [l[, these integrals can not be continued to meromorphic 
functions in the complex domain. 



V. CONCLUSION 



Thus, we have presented an algorithm for analysis of integrability of dynamical systems via the Ziglin's method 
using basically the properties of the monodromy group. As the examples show, it permits to extend the range of 
applicability of the method mainly because it resolves the problem of finding an explicit particular solution of the 
system. An important feature of the algorithm is that the trajectory obtained numerically, which is finally used for 
conclusion, is rather short and therefore can be computed with any given precision. It means that the algorithm 
indeed provides a rigorous method of computer assisted proof of non-integrability. 

We have considered some mechanical examples that were not studied before, but the algorithm also has purely 
mathematical value. Namely, it is well adapted to computation of the monodromy group (or at least some subgroup 
of it) and also represents a serious step to analysis of the differential Galois group. 

To conclude, let us mention that one of the motivations for developing constructive numerical methods of studying 
integrability for us is the qualitative analysis of the dynamical systems with delay or self-control appearing naturally 
in celestial mechanics and biological modeling. 
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